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q , Abstract 

Numerical experiments involving the interaction of two solitary waves of the 
ion acoustic plasma equations are described. An exact 2-soliton solution of the 
relevant KdV equation was fitted to the initial data, and good agreement was 
maintained throughout the entire interaction. The data demonstrates that the 
soliton interactions are virtually elastic. 
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1 Introduction 

The ion acoustic plasma equations are 

TH + (nv) x = 0, v t + I — + (f J =0, <p xx - e v + n = 0, (1) 

where <p, n, and v are respectively the electric potential, ion density and ion 
velocity. 
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We wrote a numerical code for these equations and carried out an ex- 
periment to see the interaction of two solitary waves. Solitary waves of the 
plasma equations were obtained numerically. Two solitary waves of differ- 
ent amplitudes were superposed, at a distance sufficiently far apart, and the 
evolution of the equations with this initial data was computed numerically. 
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Figure 1: Results of a numerical computation showing the interaction of two 
solitary waves for the ion acoustic plasma equations, at low amplitudes, by a 
pseudo-spectral method. The time step was dt = .008, N = 2 13 Fourier modes. 
The computation is done in a reference frame moving with speed 1.07. The speed 
of the slower wave is 1.05; while that of the larger wave is 1.1. The solid line is 
the numerical data, the dashed line an exact KdV 2 soliton solution. 

We fit the numerical data with a 2 soliton solution of the Korteweg- 
deVries (KdV) equation. Good agreement between the numerical results and 
the exact two soliton solution is maintained throughout the entire interaction. 
While numerical studies of soliton interactions in the literature, 101, 0, |Bfl, 



10], [21|, generally stress the inelastic nature of the interaction, the figures 



below show that nonelastic effects are negligible in the case of the plasma 
equations. 
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Figure 2: Upper: Trailing dispersive wave at time t = 2600; the magnitude of 
the wave is of order 10 -5 , compared with an amplitude of .15 of the slower wave. 
Bottom: If the larger wave is returned to its initial position and superposed on the 
initial wave form, the results are indisinguishable by graphical comparison; hence 
we have simply juxtaposed them here. The H 1 norm of the difference between the 
initial waveform and the translation of the final one is 6.0546 x 10~ 4 . 

In §^| we review briefly the singular perturbation argument leading to 
the KdV approximation, emphasizing the fundamental role of the Galilean 
and scaling groups. The method for fitting a two soliton solution of the 
KdV equation to the numerical data using the theoretical formulae for the 
scattering shifts of the KdV solitons is described in §||. As in the Euler 
equations, there is a wave of maximum speed and amplitude. A comparison of 
the solitary waves of the KdV equation with that of the full plasma equations 



is given in §[|. The numerical methods are discussed in §[)} Finally, in §|G], 
the results of our numerical studies will be compared with the studies in the 
literature cited above. 



2 The KdV Approximation 

The Korteweg-deVries (KdV) equation, arises as a formal singular perturba- 
tion of a number of nonlinear dispersive wave equations, such as the Euler 
equations for an inviscid fluid, and the ion acoustic equations of plasmas 



19| . Kruskal [] 1] argues that the KdV equation is the unique asymptotically 
correct model for the Euler or plasma equations under these conditions. 

Craig || has given a rigorous mathematical proof of the validity of the 
KdV and Boussinesq approximations to the Euler equations locally in time. 
His results hold for a large class of initial data, but it is unclear from his 
analysis whether the validity of the approximation is long enough to cover 
the interaction of solitary waves. Apart from this result, we know of no other 
rigorous results establishing the validity of the KdV approximation. 

The KdV equation is a solvable model. One of its special solutions is the 
so-called "multi-soliton" solution, 



d 2 / e -(e,+e k )s 
u(z,*) = 12— logdet [8 jk + ■ , (2) 

6j = u>j(x — aj — 4:Ujt), < <jj\ < ■ ■ ■ < uj n , otj e R. 

As is well known, the the solitons interact is elastically: after the interac- 
tion, they regain their original shape, and the only evidence of the interaction 
is a scattering shift. No trailing dispersive waves are generated. 

The Euler equations themselves are complicated, being a free bound- 
ary problem. The plasma equations are considerably simpler to deal with, 
both numerically and analytically, and are also an interesting physical model; 
hence, in order to investigate the validity of the KdV approximation, it makes 
sense to consider the plasma equations themselves. 

The plasma equations are Galilean invariant. This means that the equa- 
tions are the same in any Galilean frame, and we can shift to a moving frame 
of reference simply by subtracting the speed of the moving frame from the 
velocity v. In the moving frame, the velocity v' tends to — c at infinity. This 



amounts to expanding about the quiescent state n — 1, v = — c, ip = 0, where 
c is the velocity of the reference frame. 

In particular, v — — 1 at infinity in a Galilean frame moving with speed 
1, and the dispersion relation in this reference frame is found to be 

in the long wave approximation. 

The natural scaling associated with this long wave approximation is x' = 
ex, t' = e 3 t. Introducing this scaling into (0), we obtain (after division by e) 

( v 2 \ 

e 2 n t ' + (nv) x i = 0, e 2 v t > + I — + y 9 ) =0, — e 2 tp x / + e^ = n. 

This perturbation scheme is singular, since the character of the equations 
is changed when e = 0. Since only e 2 appears in these equations, we formally 
expand all quantities in powers of e 2 about n = l,v = —l,<p = 0. When 
we do this, substitute the expansions into the above equations, and collect 
terms, we find at first order that n% — V\ — (p\. 

At next order we obtain 



wi,f + {niv x ) x > + (v 2 - n 2 ) x > = 0, vi t1 ? + ( — - v 2 + <p 2 ) = 0, 

1 2 
-Vl,x'x' +<-P2+ gV 9 ! = n 2- 

The second order quantities n 2 , v 2 , and (p 2 may be eliminated from this 
system; and, dropping the primes, we obtain the Korteweg-de Vries equation 
for v\\ 

1 
Vl,t + V!Vi iX + -v 1)XXX = 0. (3) 

3 Comparison with KdV 

In order to compare the numerical data with solutions of the KdV equation, 
it is important to point out that one does not simply solve the KdV equation 
with the initial data from the plasma equations. This would produce very 



bad results, for two reasons. First, the plasma solitary wave deviates from 
the KdV solitary wave, as we show in §f|; so that if this data were taken 
as initial data for the KdV equation, secondary wave trains would develop, 
and non-elastic behavior would be observed. Moreover, there is a shift in the 



amplitudes and phases due to the soliton interaction, pl| . We now discuss 
the method of fitting an exact two soliton solution of the Korteweg deVries 
equation to the numerical data. 

First note that there is a factor of 1/2 in the KdV approximation ([3D. 
This is easily accounted for by a simple rescaling: 

v(x,t) = -u (a:, -t 

where u satisfies the KdV equation u t + uu x + u xxx = 0. 

By expanding the determinant in (|2]), and renaming the phase shifts a\, 
a 2 , the two-soliton solution of can be written 

v = §f- logr(^, 2 ), r = 1 + e~ 2ei + e" 2 * 2 + e -m+e 2+a )^ 
ax 2 

where 

9j = U)j(x — aj — (c + 2u; 2 )£), j = 1,2; a = log . 

U2 — U>i 

The two soliton solutions form a four parameter family, a%, a 2 , u>i, o; 2 . The 
linearized KdV equation at the 2 soliton solution therefore has a four dimen- 
sional null space, obtained by differentiating the equation with respect to the 
four parameters. In the perturbation series, the four parameters a%, a 2 , tO\, oji 
must be adjusted in order to eliminate the resonance terms which lie in the 



null space of the linearized operator, as discussed by Zou and Su pl| , who 
calculate the shift in these four parameters at second and third order and 
compare them with the numerical data. 

In the present case, the shifts in a\ 1 8u)2 may be determined from the 
numerical data as follows. The locations of the large and small waves before 
and after the interaction lead to four equations in four unknowns. Moreover, 
the waves are sufficiently separated before and after the interaction that we 
can apply the known formulae for the phase shifts incurred in the interaction 

The two soliton solution has the asymptotic behavior (cf. |7|) 
dx 2 



u = 6-p-^- logr ~ 6ci; 2 sech 2 (6 l i + a) + Gc^sech 2 ^, — > oo; 



and from the knowledge of the scattering shifts we find that 

u ~ 6u; 2 sech 2 1 + 6w 2 sech 2 (0 2 + ex), t — > — oo. 
Therefore, our matching conditions give 

t = : 0i = 0, 2 + a = 0; 

t = T : 0! + a = 0, 2 = 0. 

These four conditions lead to the equations 



2vr 1 , ^2 + ^1 



«i = x x = a;J~ — (c + 2u> 1 )T H lo 



u;i cj 2 — u>i 

■>, ,. _ 1 , W 2 +^l 



ot2 = x t — {c + 2u 1 )T = x 2 H log 

W 2 ^2 - Wl 

Here x ■ denote the locations of the j th wave at times t — and t = T, the 
total elapsed time, and c is the relative speed of the computation frame to 
the Galilean frame. 

The phase constants can be eliminated from these equations, and we 
obtain two equations in U\ and cj 2 : 

2cu 2 T + (-1V — log U2 + - = Sxj - cT J = 1,2, 

J UJj UJ 2 — 0J\ 

where Sxj is the total distance traversed by the j th wave. 

We solved the equations for u>i, uo 2 iteratively. As an initial guess we 
took the values obtained by matching the speeds of the two KdV waves ex- 
actly with the speeds of the soliton waves. The speed of the solitary wave 
6a; 2 sech 2 (a;(x — 2u 2 i) is 2u 2 . The speeds of the two solitary plasma waves 
(relative to the Galilean frame) are .05 and .1. Therefore, as a first approxi- 
mation, we take ui\ = a/.05/2 = .1581 and uj 2 = a/. 1/2 = .2236. 

Our data are 

T = 3800; c = -.07; x^ = 147.7349; x 2 = 50.0468; 

S Xl = -85.2439; 5x 2 = 120.4166 

We obtained ooi = .1589 and u 2 = .2231. The relative phases a,\ and a 2 are 
then determined by any of the four equations above. 
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4 Solitary waves 

The plasma equations model two-directional propagation of ion acoustic 
waves, hence the system possesses solitary waves travelling in either direction. 



Working in a frame moving with the wave, we obtain [17 



Vc 2 - 2<p (4) 



v 

n = , v = c 

c — v 

V?" = /(¥>), /(*>) = e*> - -^== (5) 



A straightforward analysis of this second order equation yields the fol- 
lowing facts, which we summarize here. Equation ([5]) has no solitary wave 
solution when < c < 1 since / is negative in a neighborhood of the origin; 
however, the origin is a stable center, and the plasma equations have a two 
parameter family (amplitude and phase) of periodic waves for all < c < 1 . 

When c > 1, the origin is a saddle point. Since f'(<p) < for all <p < 0, 
there are no homoclinic orbits in the left half phase plane. In the following 
theorem, we give a necessary and sufficient condition for the existence of a 
homoclinic orbit, which physically is a solitary wave of elevation. 

Theorem 4.1 Let c ( ~ 1.5852 ) be the positive root of the equation 

e c2 / 2 = 1 + c 2 . 

For each c G (l,c), there is a unique solitary wave of elevation to equation 
(^) which is even and real- analytic. When c = c, (|) also has a continuous 
solitary wave solution with unbounded second derivative at the origin, which 
is a weak solution of (^) in the sense of distributions. 

We omit the proof, which is not difficult. 





Figure 3. Upper: The relationship of amplitude to speed for the n and ip compo- 
nents of the plasma wave compared with that for the KdV solitary wave, which is 
linear. The amplitude of n blows up as c — > c, the maximum wave speed. Lower: 
Comparison of the single solitary KdV solution 6w 2 sech 2 u;x, (dashed), to the <p 
wave of the plasma equations (solid) . 

Theorem 4.2 For each constant c such that 1 < c < c, there is a unique 
solitary wave solution of (|l|) the form 

{1 + n(x — ct), v{x — ct), ip(x — ct)} 



where n, v, ip tend to zero at infinity. When c = c the solution is a weak 
solution in the sense of distributions. 

Left traveling waves are obtained by the reflection symmetry c i— > — c, 



v h-> —v. 



We can always make a Galilean transformation so that the asymptotic 
velocity at infinity is zero. A direct consequence of Theorems [TT] and [12] is 
the convergence of periodic orbits to solitary wave solutions. 

Corollary 4.3 Each solitary wave solution of the plasma equations is a limit 
of a sequence of periodic waves. 

5 The numerical scheme 

In this section we describe briefly the implicit pseudo-spectral scheme []20] we 
used for numerically solving the ion plasma equations. We write the plasma 
equations in a reference frame moving with speed c; we also replace n by 
n + 1, so that all three variables tend to zero at infinity. Note that this is 
not a Galilean transformation, since the velocity v still vanishes at infinity, 
and the form of the equations is changed. We then get 

n t - cn x + v x + (nv) x =0, 
fv 2 

V t -CV x + I — + if 

e v - 1 - <p + (1 - D 2 )ip =n D = d/dx. 

We integrate the first two equations with respect to time over the interval 

(ti, t 2 ), where t 2 = ti + 5t. Denote n\ = n{t\) and n 2 = n(t 2 ). We get, in the 
first equation, 

n 2 — ni+ I D(v — en + nv) dt = 0. 
Jti 

We approximate the integral using the trapezoid rule, obtaining 

1 1 

n 2 -ni + -5tD(v 2 - cn 2 + v t - aii) + -5tD(n 2 v 2 + nxVi) = 0. 

The equation in v is treated in the same way, and we obtain there 

c St I v v \ 

v 2 -v 1 - -StD( Vl +v 2 ) + jD(^- + l Pi + ^ + lp 2 \=0. 
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We combine these equations together with the equation for n 2 and (p 2 into a 
3x3 nonlinear system, which we write in the form 



(A B 0\ 


(n 2 \ 




A B 


V2 


= 


V-i o c) 


\py 





(A 1 -B \ 



A x 
\0 



-B 

/ 



W 



I — B(riiVi 



n 2 v 2 



)\ 



-Bl^ + f 



\ KM + K(<p 2 ) J 



where 

A 

and 



°^D, A 1 = 1 + C -^D, B = ^D, C = l-& 

2 2 2 



K(<p) = ±(l + <p-er). 

We found by experimentation that averaging the nonlinear term in the third 
equation produced slightly more accurate results when integrating the soli- 
tary wave. 

Inverting the 3x3 matrix on the left we obtain a nonlinear implicit 
scheme for the time step. This is solved iteratively. The size of the time 
step is constrained by the requirement that the iteration scheme converge 
sufficiently rapidly. 

The equations are integrated in a frame moving with speed 1.07. The 
speed of the larger wave is 1.1; that of the smaller wave is 1.05. The speed of 
the maximum wave is approximately 1.5852. In the moving frame, the smaller 
wave drops back, while the larger wave advances. This has the advantage of 
keeping both waves within a fixed interval throughout the interaction. 

The numerical computations are complicated by small 'discontinuities' at 
the endpoints of the interval. Though the solitary wave decays exponentially 
fast, and is of the order of 10~ 6 at the endpoints of the interval, there is nev- 
ertheless a small jump at the endpoints, if one uses the homoclinic solution, 
rather than a periodic wave nearby. Therefore, rather than computing the 
solitary wave, which is a homoclinic orbit, we instead computed a periodic 
wave very close to the homoclinic orbit. 

Numerical error introduces high frequency noise, which is transferred to 
higher modes by the nonlinear terms. For the KdV equation this effect causes 
no problems in the numerical computations, since the dispersion relation of 
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the KdV equation grows like k 3 , and acts like a high frequency cut-off in the 
implicit pseudo-spectral scheme. 

But the plasma equations have much weaker dispersion, and high fre- 
quency energies are not attenuated, with the result that the the calculations 
are corrupted over time. Because of the scaling involved, the code must be 
integrated over a very long time scale in order to see the interaction of the 
two solitary waves, in this case, to T=2600. 

To compensate for this problem, we periodically filtered the data using 
the 'sharpened raised cosine' filter (cf. [§J, p. 248): 

a{9) = ct^(35 - 84<x + 70^ - 20o$), a = -(1 + cos 9). 

The initial data was filtered; and the solution was filtered at time intervals 
of 50 units. 

The filtering, of course, removes energy from the system; but we calcu- 
lated the energy 



£ = / n(ip + cv + -v 



2^ „«, <P* 



2 



over the course of the interaction. It deviated from the original values by , 
showing that the filtering is relatively mild. 

We tested the numerical scheme by numerically intergrating a single soli- 
tary wave in a reference frame moving with the wave. The calculations 
indicate the solitary wave to be fixed point of the numerical scheme. 

6 Conclusions 

We draw two main conclusions from our study: 1) the solitary wave interac- 
tions are virtually elastic; 2) the KdV equation itself is qualitatively accurate 
in modelling the elastic interaction of solitary waves. 

Numerical computations on the second and third order terms in the sin- 
gular perturbation expansions for the Euler equations were carried out by 
Zou and Su [fJIJ]. They find that the collision is elastic up to second order; 
but the computation of the third order term shows "a secondary wave train 
trailing behind the smaller wave after the collision." They state that 

"The third-order solution contains a secondary wave train. This result 
may be considered to be consistent with the suggestion given by Benjamin 
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and Olver, (]!]]) and later Olver (||), that perfect interactions of solitary- 
waves akin to those of KdV solitons are unlikely for the complete water wave 
problem. . . " 

There is indeed a dispersive tail behind the smaller wave after the inter- 
action, as indicated in Figure 2; but it is of extremely small amplitude: ap- 
proximately 10~ 5 , compared with the amplitude of the smaller wave, which 
is approximately .15. Hence the relative amplitude of the trailing wave is 
roughly .6667 x 10~ 4 . The small parameter in our theory is u ~ .2; hence 
third order terms are of order u 6 = .64 x 10~ 4 . This is consistent with the 
dispersive wave at third order found by Zou and Su. 

A numerical study of the first and second order terms in the KdV ap- 
proximation for the ion acoustic plasma equations was carried out by Konno 
et.al. [[UJ. They found secondary waves at second order; but they did not 
cast out the secular terms in their expansion, as did Zou and Su. 

An alternative small amplitude model, for the Euler equations has been 
proposed by Benjamin et.al. ||, and is now known as the regularized long 
wave, or BBM equation. Numerical experiments by Bona, Pritchard and 
Scott H on the two soliton interaction of the BBM equation show a somewhat 
larger dispersive tail of relative magnitude 10~ 2 . Their results are confirmed 
by Kodama ||, who also shows that inelastic interactions can be expected 
in dispersive systems which are not completely integrable. 

A numerical investigation of the two soliton interaction of the full Euler 
equations was carried out by Fenton and Rienecker M. They found good 
agreement with the KdV approximation, although 'some deviations from the 
theoretical predictions were observed-the overtaking high wave grew signifi- 
cantly at the expense of the low wave, and the predicted phase shift was found 
to be only roughly described by the theory.' Comparing their results to the 
BBM equation, they state (p. 427) that their computations 'showed no trace 
of trailing waves, either from plots of the free surface or, more convincingly, 
by second differences of the point values of surface elevation'. 

Further investigation of the soliton interactions in the Euler equations 
seems warranted. Our results for the plasma equations show an excellent fit 
to the theoretical predictions of the KdV equation. Whether this is due to 
fundamental differences between the Euler and plasma equations, or due to 
the fact that the Euler equations are more difficult to integrate, is not clear. 
The Euler equations require an accurate approximation to the Dirichlet to 
Neumann map on the free surface, an issue which was not explicitly discussed 
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Ill 10. 

In our view, the inelasticity of the interaction has received undue em- 
phasis in the literature. This is misleading, since it would suggest that the 
KdV equation is not a valid model. The inelastic nature of the soliton inter- 
actions at low amplitudes are negligible, and the KdV equation models the 
interaction of overtaking solitary waves with good accuracy. 

We calculated the relative error of the KdV approximation to the v wave 
of the plasma data in the H 1 norm. We calculated the H 1 norm of the differ- 
ence between the KdV solution and the v wave, and divided by |H|#i. The 
ratio ranged from approximately .15 at t = to a maximum of .2 during the 
interaction, then dropped back to a range of .15. This is relatively high, indi- 
cating that we are operating somewhat away from the KdV approximation. 
This is clear from Figure 3, which shows a comparison of the KdV solitary 
wave with v. 

The qualitative prediction of the KdV model of elastic interactions of 
solitary waves is substantially correct. The discrepancies between the KdV 
equation and the full equations show that higher order terms in the singu- 
lar perturbation expansion must be taken into account in order to obtain a 
quantitative agreement. If anything, these two facts indicate that the elas- 
ticity of the interaction is robust, in that the interaction is virtually elastic 
even though there is a relatively large error in the KdV approximation itself. 

One natural question, considered by Moser and Sachs, [16|] is the follow- 



ing: Do the multi-soliton solutions of the KdV equation extend to solutions of 
the Euler equations on < t < oo? A positive answer to this question would 
demonstrate the global existence of a non-trivial, time dependent solution of 
the Euler equations, something which has not yet been done. 

Moser's conjecture would be established by proving the convergence of 
the singular perturbation approximation on < t < oo. This would also 1) 
establish the validity of the KdV approximation on an infinite time scale; 2) 
establish uniform estimate in time on the error term. 



Sachs |16| obtained a representation of the propagator for the linearized 
KdV equation. Using his representation, Haragus and Sattinger j?| derived 
uniform estimates for the propagator of the linearized KdV equations about 
an n-soliton solution on M x R + . The linearized equations have a 2n dimen- 
sional null space; and a Fredholm alternative for the time dependent operator 
in an appropriate function space was proved. 

In the case of the Euler and plasma equations, there is no general existence 
theorem, nor an H l estimate, as in the case of Pego and Weinstein's proof 
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13 1 of the stability of the solitary waves of the generalized KdV equation; 
hence their analysis does not directly apply here. We anticipate that a Nash- 
Moser type of implicit function theorem [[OJ, [[18] may be needed. A natural 
space to work in is the space of functions analytic in a strip in the complex 
x plane. Estimates for the linearized operator in this class were obtained in 
which were uniform on < t < oo. The implications of this are that 
the perturbation series itself does not converge but is only asymptotically 
valid. The actual solution would be C°°, and accurate to all orders, but not 
analytic in e. 
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